# Replication code for 'Policy Preferences in a Post-War Environment'

# NOTE: when replicating the graphical analyses, ensure that the graphing/plotting window is sufficiently large.

# clear working directory
rm(list=ls()) 

# set the appropriate working directory; Code.R and Data.csv should be located in the same directory 
# setwd(...) 

# install all required packages
# install.packages(...)

# load all required packages
library(interplot)
library(ordinal)
library(memisc)
library(nnet)
library(coefplot)
library(effects)
library(pmr)
library(gridExtra)
library(survey)
library(AER)

# load data 
data <- read.csv("Data.csv")

# set reference level for Employment_Status and order the ordinal variables
data$Employment_Status <- relevel(data$Employment_Status, ref="Not in labor force")
data$Child_Care_OP <- factor(data$Child_Care_OP, levels=c("Least Important","Second Least Important","Second Most Important","Most Important"),ordered=TRUE)
data$Education_OP <- factor(data$Education_OP, levels=c("Least Important","Second Least Important","Second Most Important","Most Important"),ordered=TRUE)
data$Infrastructure_OP <- factor(data$Infrastructure_OP, levels=c("Least Important","Second Least Important","Second Most Important","Most Important"),ordered=TRUE)
data$Policing_OP <- factor(data$Policing_OP, levels=c("Least Important","Second Least Important","Second Most Important","Most Important"),ordered=TRUE)

# create datasets which will be used in the analyses
data_C_VP <- data[data$Condition != "Identity_Prime",]
data_C_IP <- data[data$Condition != "Violence_Prime",]
data_IP_VP <- data[data$Condition != "Control",]
data_C_VP$Treatment <- ifelse(data_C_VP$Condition == "Violence_Prime", 1, 0)
data_C_IP$Treatment <- ifelse(data_C_IP$Condition == "Identity_Prime", 1, 0)
data_IP_VP$Treatment <- ifelse(data_IP_VP$Condition == "Violence_Prime", 1, 0)
data_Control <- data[data$Condition == "Control",]
data_IP <- data[data$Condition == "Identity_Prime",]
data_VP <- data[data$Condition == "Violence_Prime",]

# Policy Preferences in a Post-War Environment
# Main text

# Figure 1: Safety Threat and Policy Preferences
# Left Panel
# Baseline: Control
# Treatment: Violence Prime
OLS_Left <- lm(Policing_OLS ~ Treatment, data=data_C_VP)
Oprobit_Left <- clm(Policing_OP ~ Treatment, data=data_C_VP, Hess = TRUE, link = "probit")
Probit_Left <- glm(Policing_Probit ~ Treatment, data=data_C_VP, family=binomial(link="probit"))
par(mfrow=c(1,3), mar=c(5.1,9.5,4.1,0)) 
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
axis(2, at=c(0.25, 0.5, 0.75), lwd=1, 
     labels = c("Probit", 
                "Ordered Probit",
                "OLS"), las=TRUE, cex.axis=1.5)
summary(OLS_Left)
summary(Oprobit_Left)
summary(Probit_Left)
points(x=-0.05168, y=0.75, pch=19, lwd=2)
points(x=-0.04705, y=0.5, pch=19, lwd=2)
points(x=-0.07759, y=0.25, pch=19, lwd=2)
confint(OLS_Left)
confint(Oprobit_Left)
confint(Probit_Left)
segments(x0= -0.2098206, 0.75, x1 = 0.1064614, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2019169, 0.5, x1 = 0.1077979, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2670765, 0.25, x1 = 0.1118301, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

# Middle Panel
# Baseline: Control
# Treatment: Identity Prime
OLS_Middle <- lm(Policing_OLS ~ Treatment, data=data_C_IP)
Oprobit_Middle <- clm(Policing_OP ~ Treatment, data=data_C_IP, Hess = TRUE, link = "probit")
Probit_Middle <- glm(Policing_Probit ~ Treatment, data=data_C_IP, family=binomial(link="probit"))
par(mar=c(5.1,5,4.1,4.5))
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_Middle)
summary(Oprobit_Middle)
summary(Probit_Middle)
points(x=0.03175, y=0.75, pch=19, lwd=2)
points(x=0.02888, y=0.5, pch=19, lwd=2)
points(x=-0.05775, y=0.25, pch=19, lwd=2)
confint(OLS_Middle)
confint(Oprobit_Middle)
confint(Probit_Middle)
segments(x0= -0.125044, 0.75, x1 = 0.188536, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.1267198, 0.5, x1 = 0.1844832, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2479393, 0.25, x1 = 0.1323955, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Identity Prime"), line=1, cex.main=1.5)

# Right Panel
# Baseline: Identity Prime
# Treatment: Violence Prime
OLS_Right <- lm(Policing_OLS ~ Treatment, data=data_IP_VP)
Oprobit_Right <- clm(Policing_OP ~ Treatment, data=data_IP_VP, Hess = TRUE, link = "probit")
Probit_Right <- glm(Policing_Probit ~ Treatment, data=data_IP_VP, family=binomial(link="probit"))
par(mar=c(5.1,0.5,4.1,9))
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_Right)
summary(Oprobit_Right)
summary(Probit_Right)
points(x=-0.08343, y=0.75, pch=19, lwd=2)
points(x=-0.07910, y=0.5, pch=19, lwd=2)
points(x=-0.01985, y=0.25, pch=19, lwd=2)
confint(OLS_Right)
confint(Oprobit_Right)
confint(Probit_Right)
segments(x0= -0.2358173, 0.75, x1 = 0.06896609, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2316126, 0.5, x1 = 0.07339672, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2080324, 0.25, x1 = 0.1683300, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Identity Prime"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

# Policy Preferences in a Post-War Environment
# Supplementary Appendix (SA)

# SA3: Summary Statistics
# Table SA3.1: Summary Statistics
# number of respondents
dim(data)[1] 

# number of respondents assigned to each condition
table(data$Condition) 
# percentage of respondents assigned to each condition
(table(data$Condition)/1125)*100 
# no non-response for condition assignment
sum(!complete.cases(data$Condition)) 

# number of respondents from each ethnic group
table(data$Ethnicity) 
# percentage of respondents from each ethnic group
(table(data$Ethnicity)/1125)*100 
# no non-response for ethnicity
sum(!complete.cases(data$Ethnicity)) 

# 0 = number of women, 1 = number of men
table(data$Male) 
# 0 = percentage of women, 1 = percentage of men
(table(data$Male)/1125)*100 
# no non-response for gender
sum(!complete.cases(data$Male)) 

# 0 = number unmarried, 1 = number married
table(data$Married) 
# 0 = percentage unmarried, 1 = percentage married
(table(data$Married)/1125)*100 
# non-response for marital status (number)
sum(!complete.cases(data$Married)) 
# non-response for marital status (percentage)
(sum(!complete.cases(data$Married))/1125)*100 

# number unemployed, employed, and not in labor force
table(data$Employment_Status) 
# percentage unemployed, employed, and not in labor force
(table(data$Employment_Status)/1125)*100 
# non-response for employment status (number)
sum(!complete.cases(data$Employment_Status)) 
# non-response for employment status (percentage)
(sum(!complete.cases(data$Employment_Status))/1125)*100 

# mean age
mean(data$Age) 
# median age
median(data$Age) 
# no non-response for age
sum(!complete.cases(data$Age)) 

# mean education
mean(data$Education) 
# median education
median(data$Education) 
# no non-response for education
sum(!complete.cases(data$Education)) 

# SA5: Balance Tests
# Table SA5.1: Balance Tests on Respondent Characteristics
# Model 1: Violence Prime
summary(modSA5.1.1 <- glm(Treatment ~ Ethnicity + Male + Age + Education, data=data_C_VP, family=binomial(link = "probit")))
# extract the number of observations
nobs(modSA5.1.1)

# Model 2: Violence Prime
summary(modSA5.1.2 <- glm(Treatment ~ Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_VP, family=binomial(link = "probit")))
# extract the number of observations
nobs(modSA5.1.2)

# Model 3: Identity Prime
summary(modSA5.1.3 <- glm(Treatment ~ Ethnicity + Male + Age + Education, data=data_C_IP, family=binomial(link = "probit")))
# extract the number of observations
nobs(modSA5.1.3)

# Model 4: Identity Prime
summary(modSA5.1.4 <- glm(Treatment ~ Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_IP, family=binomial(link = "probit")))
# extract the number of observations
nobs(modSA5.1.4)

# Model 5: Violence Prime
summary(modSA5.1.5 <- glm(Treatment ~ Ethnicity + Male + Age + Education , data=data_IP_VP, family=binomial(link = "probit")))
# extract the number of observations
nobs(modSA5.1.5)

# Model 6: Violence Prime
summary(modSA5.1.6 <- glm(Treatment ~ Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_IP_VP, family=binomial(link = "probit")))
# extract the number of observations
nobs(modSA5.1.6)

# SA7: Distribution of Respnses Across Experimental Conditions
# Table SA7.1: Distribution of Responses across Full Sample
# number of respondents, child care
table(data$Child_Care_OP)
# percentage of respondents, child care
(table(data$Child_Care_OP)/sum(table(data$Child_Care_OP)))*100

# number of respondents, education
table(data$Education_OP)
# percentage of respondents, education
(table(data$Education_OP)/sum(table(data$Education_OP)))*100

# number of respondents, infrastructure
table(data$Infrastructure_OP)
# percentage of respondents, infrastructure
(table(data$Infrastructure_OP)/sum(table(data$Infrastructure_OP)))*100

# number of respondents, policing
table(data$Policing_OP)
# percentage of respondents, policing
(table(data$Policing_OP)/sum(table(data$Policing_OP)))*100

# Table SA7.2: Distribution of Responses across Control Condition
# number of respondents, child care
table(data_Control$Child_Care_OP)
# percentage of respondents, child care
(table(data_Control$Child_Care_OP)/sum(table(data_Control$Child_Care_OP)))*100

# number of respondents, education
table(data_Control$Education_OP)
# percentage of respondents, education
(table(data_Control$Education_OP)/sum(table(data_Control$Education_OP)))*100

# number of respondents, infrastructure
table(data_Control$Infrastructure_OP)
# percentage of respondents, infrastructure
(table(data_Control$Infrastructure_OP)/sum(table(data_Control$Infrastructure_OP)))*100

# number of respondents, policing
table(data_Control$Policing_OP)
# percentage of respondents, policing
(table(data_Control$Policing_OP)/sum(table(data_Control$Policing_OP)))*100

# Table SA7.3: Distribution of Responses across Violence Prime Condition
# number of respondents, child care
table(data_VP$Child_Care_OP)
# percentage of respondents, child care
(table(data_VP$Child_Care_OP)/sum(table(data_VP$Child_Care_OP)))*100

# number of respondents, education
table(data_VP$Education_OP)
# percentage of respondents, education
(table(data_VP$Education_OP)/sum(table(data_VP$Education_OP)))*100

# number of respondents, infrastructure
table(data_VP$Infrastructure_OP)
# percentage of respondents, infrastructure
(table(data_VP$Infrastructure_OP)/sum(table(data_VP$Infrastructure_OP)))*100

# number of respondents, policing
table(data_VP$Policing_OP)
# percentage of respondents, policing
(table(data_VP$Policing_OP)/sum(table(data_VP$Policing_OP)))*100

# Table SA7.4: Distribution of Responses across Identity Prime Condition
# number of respondents, child care
table(data_IP$Child_Care_OP)
# percentage of respondents, child care
(table(data_IP$Child_Care_OP)/sum(table(data_IP$Child_Care_OP)))*100

# number of respondents, education
table(data_IP$Education_OP)
# percentage of respondents, education
(table(data_IP$Education_OP)/sum(table(data_IP$Education_OP)))*100

# number of respondents, infrastructure
table(data_IP$Infrastructure_OP)
# percentage of respondents, infrastructure
(table(data_IP$Infrastructure_OP)/sum(table(data_IP$Infrastructure_OP)))*100

# number of respondents, policing
table(data_IP$Policing_OP)
# percentage of respondents, policing
(table(data_IP$Policing_OP)/sum(table(data_IP$Policing_OP)))*100

# SA8: Treatment Effects Across Moderating Variables
# Figure SA8.1: Safety Threat and Policy Preferences across Ethnic Groups
# Bosniaks
par(mfrow=c(1,3), mar=c(5.1,9.5,4.1,0)) 
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
axis(2, at=c(0.25, 0.5, 0.75), lwd=1, 
     labels = c("Probit", 
                "Ordered Probit",
                "OLS"), las=TRUE, cex.axis=1.5)
summary(OLS_C_VP_B <- lm(Policing_OLS ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Bosniak",])) 
summary(Oprobit_C_VP_B <- clm(Policing_OP ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Bosniak",], Hess = TRUE, link = "probit")) 
summary(Probit_C_VP_B <- glm(Policing_Probit ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Bosniak",], family=binomial(link="probit"))) 
points(x=-0.02201, y=0.75, pch=19, lwd=2)
points(x=-0.02554, y=0.5, pch=19, lwd=2)
points(x=-0.06965, y=0.25, pch=19, lwd=2)
confint(OLS_C_VP_B)
confint(Oprobit_C_VP_B)
confint(Probit_C_VP_B)
segments(x0= -0.2335907, 0.75, x1 = 0.1895618, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2384787, 0.5, x1 = 0.1873474, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.3239726, 0.25, x1 = 0.1845646, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

par(mar=c(5.1,5,4.1,4.5))
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_C_IP_B <- lm(Policing_OLS ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Bosniak",])) 
summary(Oprobit_C_IP_B <- clm(Policing_OP ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Bosniak",], Hess = TRUE, link = "probit")) 
summary(Probit_C_IP_B <- glm(Policing_Probit ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Bosniak",], family=binomial(link="probit"))) 
points(x=0.002492, y=0.75, pch=19, lwd=2)
points(x=-0.004434, y=0.5, pch=19, lwd=2)
points(x=-0.07804, y=0.25, pch=19, lwd=2)
confint(OLS_C_IP_B)
confint(Oprobit_C_IP_B)
confint(Probit_C_IP_B)
segments(x0= -0.2060701, 0.75, x1 = 0.2110541, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2160168, 0.5, x1 = 0.2070991, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.3310300, 0.25, x1 = 0.1748723, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Identity Prime"), line=1, cex.main=1.5)

par(mar=c(5.1,0.5,4.1,9))
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_IP_VP_B <- lm(Policing_OLS ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Bosniak",])) 
summary(Oprobit_IP_VP_B <- clm(Policing_OP ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Bosniak",], Hess = TRUE, link = "probit")) 
summary(Probit_IP_VP_B <- glm(Policing_Probit ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Bosniak",], family=binomial(link="probit"))) 
points(x=-0.02451, y=0.75, pch=19, lwd=2)
points(x=-0.02185, y=0.5, pch=19, lwd=2)
points(x=0.008387, y=0.25, pch=19, lwd=2)
confint(OLS_IP_VP_B)
confint(Oprobit_IP_VP_B)
confint(Probit_IP_VP_B)
segments(x0= -0.2278435, 0.75, x1 = 0.1788308, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2296311, 0.5, x1 = 0.1859261, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2420633, 0.25, x1 = 0.2588115, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Identity Prime"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

# Croats
par(mfrow=c(1,3), mar=c(5.1,9.5,4.1,0)) 
frame()
plot.window(xlim=c(-1.5, 1.5), ylim=c(0, 1))
axis(1, at=c(-1.5, -0.75, 0, 0.75, 1.5), lwd=1, cex.axis=1.5)
box(lwd = 1)
axis(2, at=c(0.25, 0.5, 0.75), lwd=1, 
     labels = c("Probit", 
                "Ordered Probit",
                "OLS"), las=TRUE, cex.axis=1.5)
summary(OLS_C_VP_C <- lm(Policing_OLS ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Croat",])) 
summary(Oprobit_C_VP_C <- clm(Policing_OP ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Croat",], Hess = TRUE, link = "probit")) 
summary(Probit_C_VP_C <- glm(Policing_Probit ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Croat",], family=binomial(link="probit"))) 
points(x=-0.2292, y=0.75, pch=19, lwd=2)
points(x=-0.1809, y=0.5, pch=19, lwd=2)
points(x=-0.5039, y=0.25, pch=19, lwd=2)
confint(OLS_C_VP_C)
confint(Oprobit_C_VP_C)
confint(Probit_C_VP_C)
segments(x0= -0.7391316, 0.75, x1 = 0.2807983, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.6410022, 0.5, x1 = 0.2787525, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -1.0995454, 0.25, x1 = 0.07476581, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

par(mar=c(5.1,5,4.1,4.5))
frame()
plot.window(xlim=c(-1.5, 1.5), ylim=c(0, 1))
axis(1, at=c(-1.5, -0.75, 0, 0.75, 1.5), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_C_IP_C <- lm(Policing_OLS ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Croat",])) 
summary(Oprobit_C_IP_C <- clm(Policing_OP ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Croat",], Hess = TRUE, link = "probit")) 
summary(Probit_C_IP_C <- glm(Policing_Probit ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Croat",], family=binomial(link="probit"))) 
points(x=0.06137, y=0.75, pch=19, lwd=2)
points(x=0.04184, y=0.5, pch=19, lwd=2)
points(x=-0.1015, y=0.25, pch=19, lwd=2)
confint(OLS_C_IP_C)
confint(Oprobit_C_IP_C)
confint(Probit_C_IP_C)
segments(x0= -0.4856288, 0.75, x1 = 0.6083766, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.4323874, 0.5, x1 = 0.516099, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.6661819, 0.25, x1 = 0.45832994, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Identity Prime"), line=1, cex.main=1.5)

par(mar=c(5.1,0.5,4.1,9))
frame()
plot.window(xlim=c(-1.5, 1.5), ylim=c(0, 1))
axis(1, at=c(-1.5, -0.75, 0, 0.75, 1.5), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_IP_VP_C <- lm(Policing_OLS ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Croat",])) 
summary(Oprobit_IP_VP_C <- clm(Policing_OP ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Croat",], Hess = TRUE, link = "probit")) 
summary(Probit_IP_VP_C <- glm(Policing_Probit ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Croat",], family=binomial(link="probit"))) 
points(x=-0.2905, y=0.75, pch=19, lwd=2)
points(x=-0.2531, y=0.5, pch=19, lwd=2)
points(x=-0.4024, y=0.25, pch=19, lwd=2)
confint(OLS_IP_VP_C)
confint(Oprobit_IP_VP_C)
confint(Probit_IP_VP_C)
segments(x0= -0.8117716, 0.75, x1 = 0.2306905, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.7400994, 0.5, x1 = 0.2335079, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -1.0335153, 0.25, x1 = 0.2176774, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Identity Prime"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

# Serbs
par(mfrow=c(1,3), mar=c(5.1,9.5,4.1,0)) 
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
axis(2, at=c(0.25, 0.5, 0.75), lwd=1, 
     labels = c("Probit", 
                "Ordered Probit",
                "OLS"), las=TRUE, cex.axis=1.5)
summary(OLS_C_VP_S <- lm(Policing_OLS ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Serb",])) 
summary(Oprobit_C_VP_S <- clm(Policing_OP ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Serb",], Hess = TRUE, link = "probit")) 
summary(Probit_C_VP_S <- glm(Policing_Probit ~ Treatment, data=data_C_VP[data_C_VP$Ethnicity == "Serb",], family=binomial(link="probit"))) 
points(x=-0.04486, y=0.75, pch=19, lwd=2)
points(x=-0.03819, y=0.5, pch=19, lwd=2)
points(x=0.05807, y=0.25, pch=19, lwd=2)
confint(OLS_C_VP_S)
confint(Oprobit_C_VP_S)
confint(Probit_C_VP_S)
segments(x0= -0.3095906, 0.75, x1 = 0.2198758, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2996314, 0.5, x1 = 0.2232226, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.275048, 0.25, x1 = 0.3927981, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

par(mar=c(5.1,5,4.1,4.5))
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_C_IP_S <- lm(Policing_OLS ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Serb",])) 
summary(Oprobit_C_IP_S <- clm(Policing_OP ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Serb",], Hess = TRUE, link = "probit")) 
summary(Probit_C_IP_S <- glm(Policing_Probit ~ Treatment, data=data_C_IP[data_C_IP$Ethnicity == "Serb",], family=binomial(link="probit"))) 
points(x=0.04997, y=0.75, pch=19, lwd=2)
points(x=0.0562, y=0.5, pch=19, lwd=2)
points(x=0.004374, y=0.25, pch=19, lwd=2)
confint(OLS_C_IP_S)
confint(Oprobit_C_IP_S)
confint(Probit_C_IP_S)
segments(x0= -0.2099504, 0.75, x1 = 0.3098927, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2092816, 0.5, x1 = 0.3217272, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.3376641, 0.25, x1 = 0.3471180, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Control"), line=2.5, cex.main=1.5)
title(expression("Treatment: Identity Prime"), line=1, cex.main=1.5)

par(mar=c(5.1,0.5,4.1,9))
frame()
plot.window(xlim=c(-0.4, 0.4), ylim=c(0, 1))
axis(1, at=c(-0.4, -0.2, 0, 0.2, 0.4), lwd=1, cex.axis=1.5)
box(lwd = 1)
summary(OLS_IP_VP_S <- lm(Policing_OLS ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Serb",])) 
summary(Oprobit_IP_VP_S <- clm(Policing_OP ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Serb",], Hess = TRUE, link = "probit")) 
summary(Probit_IP_VP_S <- glm(Policing_Probit ~ Treatment, data=data_IP_VP[data_IP_VP$Ethnicity == "Serb",], family=binomial(link="probit"))) 
points(x=-0.09483, y=0.75, pch=19, lwd=2)
points(x=-0.09701, y=0.5, pch=19, lwd=2)
points(x=0.0537, y=0.25, pch=19, lwd=2)
confint(OLS_IP_VP_S)
confint(Oprobit_IP_VP_S)
confint(Probit_IP_VP_S)
segments(x0= -0.3479476, 0.75, x1 = 0.1582906, y1 = 0.75,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.35174, 0.5, x1 = 0.1576296, y1 = 0.5,
         lty = par("lty"), lwd = par("lwd"))
segments(x0= -0.2710110, 0.25, x1 = 0.3793180, y1 = 0.25,
         lty = par("lty"), lwd = par("lwd"))
segments(x0=0, x1=0, y0=-.04, y1=1.04, lty=2)
title(expression("Baseline: Identity Prime"), line=2.5, cex.main=1.5)
title(expression("Treatment: Violence Prime"), line=1, cex.main=1.5)

# Figure SA8.2: Safety Threat and Policy Preferences across Local Violence Severity
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*Casualty, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="Casualty") + ggtitle("OLS") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*Casualty, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="Casualty") + ggtitle("Ordered Probit") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*Casualty, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="Casualty") + ggtitle("Probit") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*Casualty, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="Casualty") + ggtitle("OLS") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*Casualty, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="Casualty") + ggtitle("Ordered Probit") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*Casualty, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="Casualty") + ggtitle("Probit") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*Casualty, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="Casualty") + ggtitle("OLS") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*Casualty, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="Casualty") + ggtitle("Ordered Probit") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*Casualty, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="Casualty") + ggtitle("Probit") +
  xlab("Casualty") + ylab("Treatment Effect") + ylim(c(-1.25, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Figure SA8.3: Safety Threat and Policy Preferences across Age
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*Age, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="Age") + ggtitle("OLS") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*Age, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="Age") + ggtitle("Ordered Probit") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*Age, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="Age") + ggtitle("Probit") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*Age, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="Age") + ggtitle("OLS") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*Age, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="Age") + ggtitle("Ordered Probit") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*Age, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="Age") + ggtitle("Probit") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*Age, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="Age") + ggtitle("OLS") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*Age, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="Age") + ggtitle("Ordered Probit") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*Age, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="Age") + ggtitle("Probit") +
  xlab("Respondent's Age") + ylab("Treatment Effect") + ylim(c(-1, 1)) + xlim(c(18, 85)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Figure SA8.4: Safety Threat and Policy Preferences across Age Status
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*Adult_During_War, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="Adult_During_War") + ggtitle("OLS") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*Adult_During_War, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="Adult_During_War") + ggtitle("Ordered Probit") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*Adult_During_War, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="Adult_During_War") + ggtitle("Probit") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*Adult_During_War, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="Adult_During_War") + ggtitle("OLS") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*Adult_During_War, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="Adult_During_War") + ggtitle("Ordered Probit") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*Adult_During_War, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="Adult_During_War") + ggtitle("Probit") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*Adult_During_War, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="Adult_During_War") + ggtitle("OLS") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*Adult_During_War, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="Adult_During_War") + ggtitle("Ordered Probit") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*Adult_During_War, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="Adult_During_War") + ggtitle("Probit") +
  xlab("Respondent's Age Status During War") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") +
  scale_x_continuous(breaks = c(0,1), labels = c("Not Adult", "Adult"))
# Message indicates that we are replacing the original numeric x-axis with a new, character x-axis: "Not Adult" and "Adult"

# Figure SA8.5: Safety Threat and Policy Preferences across Standard Ethnic Fractionalization Index
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*FRAC, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="FRAC") + ggtitle("OLS") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*FRAC, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="FRAC") + ggtitle("Ordered Probit") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*FRAC, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="FRAC") + ggtitle("Probit") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*FRAC, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="FRAC") + ggtitle("OLS") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*FRAC, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="FRAC") + ggtitle("Ordered Probit") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*FRAC, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="FRAC") + ggtitle("Probit") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*FRAC, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="FRAC") + ggtitle("OLS") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*FRAC, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="FRAC") + ggtitle("Ordered Probit") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*FRAC, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="FRAC") + ggtitle("Probit") +
  xlab("Fractionalization Index") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Figure SA8.6: Safety Threat and Policy Preferences across Modified Ethnic Fractionalization Index (Vignette Groups)
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*FRAC_Respondent_Vignette, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("OLS") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*FRAC_Respondent_Vignette, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("Ordered Probit") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*FRAC_Respondent_Vignette, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("Probit") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*FRAC_Respondent_Vignette, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("OLS") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*FRAC_Respondent_Vignette, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("Ordered Probit") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*FRAC_Respondent_Vignette, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("Probit") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*FRAC_Respondent_Vignette, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("OLS") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*FRAC_Respondent_Vignette, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("Ordered Probit") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*FRAC_Respondent_Vignette, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="FRAC_Respondent_Vignette") + ggtitle("Probit") +
  xlab("FRAC: Respondent and Vignette Groups") + ylab("Treatment Effect") + ylim(c(-1.25, 1.25)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted")

# Figure SA8.7: Safety Threat and Policy Preferences across Modified Ethnic Fractionalization Index (Respondent and all Out-Group Members)
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*FRAC_Respondent_OutGroups, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("OLS") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*FRAC_Respondent_OutGroups, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("Ordered Probit") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*FRAC_Respondent_OutGroups, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("Probit") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*FRAC_Respondent_OutGroups, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("OLS") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*FRAC_Respondent_OutGroups, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("Ordered Probit") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*FRAC_Respondent_OutGroups, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("Probit") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*FRAC_Respondent_OutGroups, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("OLS") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*FRAC_Respondent_OutGroups, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("Ordered Probit") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*FRAC_Respondent_OutGroups, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="FRAC_Respondent_OutGroups") + ggtitle("Probit") +
  xlab("FRAC: Respondent and all Out-Groups") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Figure SA8.8: Safety Threat and Policy Preferences across Co-Ethnics as Percentage of Municipal Population
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP)) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP, Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP, family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1, 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Figure SA8.9: Safety Threat and Policy Preferences across Co-Ethnics as Percentage of Municipal Population, Respondents from Republika Srpska
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP[data_C_VP$RS == 1,])) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP[data_C_VP$RS == 1,], Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP[data_C_VP$RS == 1,], family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP[data_C_IP$RS == 1,])) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP[data_C_IP$RS == 1,], Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP[data_C_IP$RS == 1,], family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP[data_IP_VP$RS == 1,])) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP[data_IP_VP$RS == 1,], Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP[data_IP_VP$RS == 1,], family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Figure SA8.10: Safety Threat and Policy Preferences across Co-Ethnics as Percentage of Municipal Population, Respondents from FBiH/Brcko
# Baseline = Control, Treatment = Violence Prime
summary(OLS_C_VP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP[data_C_VP$FBiH_Brcko == 1,])) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_VP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP[data_C_VP$FBiH_Brcko == 1,], Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_VP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_C_VP[data_C_VP$FBiH_Brcko == 1,], family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Control, Treatment = Identity Prime
summary(OLS_C_IP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP[data_C_IP$FBiH_Brcko == 1,])) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_C_IP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP[data_C_IP$FBiH_Brcko == 1,], Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_C_IP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_C_IP[data_C_IP$FBiH_Brcko == 1,], family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_C_IP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# Baseline = Identity Prime, Treatment = Violence Prime
summary(OLS_IP_VP <- lm(Policing_OLS ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP[data_IP_VP$FBiH_Brcko == 1,])) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(OLS_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("OLS") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Oprobit_IP_VP <- polr(Policing_OP ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP[data_IP_VP$FBiH_Brcko == 1,], Hess = TRUE, method = "probit")) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Oprobit_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Ordered Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

summary(Probit_IP_VP <- glm(Policing_Probit ~ Treatment*Respondent_Per_CoEthnic, data=data_IP_VP[data_IP_VP$FBiH_Brcko == 1,], family=binomial(link="probit"))) 
theme_update(plot.title = element_text(hjust = 0.5))
interplot(Probit_IP_VP, var1="Treatment", var2="Respondent_Per_CoEthnic") + ggtitle("Probit") +
  xlab("% Co-Ethnic") + ylab("Treatment Effect") + ylim(c(-1.5, 1.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_hline(yintercept=c(0), linetype="dotted") 

# SA9: Treatment Effects Across Moderating Variables
# Table SA9.1: Safety Threat and Policy Preferences, OLS Models
# Model 1
summary(SA9.1.1 <- lm(Policing_OLS ~ Treatment, data=data_C_VP))
# extract the number of observations
nobs(SA9.1.1)
# extract the AIC
AIC(SA9.1.1)

# Model 2
summary(SA9.1.2 <- lm(Policing_OLS ~ Treatment, data=data_C_IP))
# extract the number of observations
nobs(SA9.1.2)
# extract the AIC
AIC(SA9.1.2)

# Model 3
summary(SA9.1.3 <- lm(Policing_OLS ~ Treatment, data=data_IP_VP))
# extract the number of observations
nobs(SA9.1.3)
# extract the AIC
AIC(SA9.1.3)

# Table SA9.2: Safety Threat and Policy Preferences, Ordered Probit Models
# Model 1
summary(SA9.2.1 <- clm(Policing_OP ~ Treatment, data=data_C_VP, link = "probit"))
# extract the number of observations
nobs(SA9.2.1)
# extract the AIC
AIC(SA9.2.1)

# Model 2
summary(SA9.2.2 <- clm(Policing_OP ~ Treatment, data=data_C_IP, link = "probit"))
# extract the number of observations
nobs(SA9.2.2)
# extract the AIC
AIC(SA9.2.2)

# Model 3
summary(SA9.2.3 <- clm(Policing_OP ~ Treatment, data=data_IP_VP, link = "probit"))
# extract the number of observations
nobs(SA9.2.3)
# extract the AIC
AIC(SA9.2.3)

# Table SA9.3: Safety Threat and Policy Preferences, Probit Models
# Model 1
summary(SA9.3.1 <- glm(Policing_Probit ~ Treatment, data=data_C_VP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.3.1)
# extract the AIC
AIC(SA9.3.1)

# Model 2
summary(SA9.3.2 <- glm(Policing_Probit ~ Treatment, data=data_C_IP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.3.2)
# extract the AIC
AIC(SA9.3.2)

# Model 3
summary(SA9.3.3 <- glm(Policing_Probit ~ Treatment, data=data_IP_VP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.3.3)
# extract the AIC
AIC(SA9.3.3)

# Table SA9.4: Safety Threat and Policy Preferences, OLS Models, Controlling for Respondent Characteristics
# Model 1
summary(SA9.4.1 <- lm(Policing_OLS ~ Treatment + Ethnicity + Male + Age + Education, data=data_C_VP))
# extract the number of observations
nobs(SA9.4.1)
# extract the AIC
AIC(SA9.4.1)

# Model 2
summary(SA9.4.2 <- lm(Policing_OLS ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_VP))
# extract the number of observations
nobs(SA9.4.2)
# extract the AIC
AIC(SA9.4.2)

# Model 3
summary(SA9.4.3 <- lm(Policing_OLS ~ Treatment + Ethnicity + Male + Age + Education, data=data_C_IP))
# extract the number of observations
nobs(SA9.4.3)
# extract the AIC
AIC(SA9.4.3)

# Model 4
summary(SA9.4.4 <- lm(Policing_OLS ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_IP))
# extract the number of observations
nobs(SA9.4.4)
# extract the AIC
AIC(SA9.4.4)

# Model 5
summary(SA9.4.5 <- lm(Policing_OLS ~ Treatment + Ethnicity + Male + Age + Education, data=data_IP_VP))
# extract the number of observations
nobs(SA9.4.5)
# extract the AIC
AIC(SA9.4.5)

# Model 6
summary(SA9.4.6 <- lm(Policing_OLS ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_IP_VP))
# extract the number of observations
nobs(SA9.4.6)
# extract the AIC
AIC(SA9.4.6)

# Table SA9.5: Safety Threat and Policy Preferences, Ordered Probit Models, Controlling for Respondent Characteristics
# Model 1
summary(SA9.5.1 <- clm(Policing_OP ~ Treatment + Ethnicity + Male + Age + Education, data=data_C_VP, link = "probit"))
# extract the number of observations
nobs(SA9.5.1)
# extract the AIC
AIC(SA9.5.1)

# Model 2
summary(SA9.5.2 <- clm(Policing_OP ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_VP, link = "probit"))
# extract the number of observations
nobs(SA9.5.2)
# extract the AIC
AIC(SA9.5.2)

# Model 3
summary(SA9.5.3 <- clm(Policing_OP ~ Treatment + Ethnicity + Male + Age + Education, data=data_C_IP, link = "probit"))
# extract the number of observations
nobs(SA9.5.3)
# extract the AIC
AIC(SA9.5.3)

# Model 4
summary(SA9.5.4 <- clm(Policing_OP ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_IP, link = "probit"))
# extract the number of observations
nobs(SA9.5.4)
# extract the AIC
AIC(SA9.5.4)

# Model 5
summary(SA9.5.5 <- clm(Policing_OP ~ Treatment + Ethnicity + Male + Age + Education, data=data_IP_VP, link = "probit"))
# extract the number of observations
nobs(SA9.5.5)
# extract the AIC
AIC(SA9.5.5)

# Model 6
summary(SA9.5.6 <- clm(Policing_OP ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_IP_VP, link = "probit"))
# extract the number of observations
nobs(SA9.5.6)
# extract the AIC
AIC(SA9.5.6)

# Table SA9.6: Safety Threat and Policy Preferences, Probit Models, Controlling for Respondent Characteristics
# Model 1
summary(SA9.6.1 <- glm(Policing_Probit ~ Treatment + Ethnicity + Male + Age + Education, data=data_C_VP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.6.1)
# extract the AIC
AIC(SA9.6.1)

# Model 2
summary(SA9.6.2 <- glm(Policing_Probit ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_VP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.6.2)
# extract the AIC
AIC(SA9.6.2)

# Model 3
summary(SA9.6.3 <- glm(Policing_Probit ~ Treatment + Ethnicity + Male + Age + Education, data=data_C_IP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.6.3)
# extract the AIC
AIC(SA9.6.3)

# Model 4
summary(SA9.6.4 <- glm(Policing_Probit ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_C_IP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.6.4)
# extract the AIC
AIC(SA9.6.4)

# Model 5
summary(SA9.6.5 <- glm(Policing_Probit ~ Treatment + Ethnicity + Male + Age + Education, data=data_IP_VP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.6.5)
# extract the AIC
AIC(SA9.6.5)

# Model 6
summary(SA9.6.6 <- glm(Policing_Probit ~ Treatment + Ethnicity + Male + Age + Education + Married + Employment_Status, data=data_IP_VP, family=binomial(link="probit")))
# extract the number of observations
nobs(SA9.6.6)
# extract the AIC
AIC(SA9.6.6)

# ANALYSIS WITH SURVEY WEIGHTS
data$Employment_Status_Census <- relevel(data$Employment_Status_Census, "Not in labor force")

data <- data[complete.cases(data$Male_Census),]
data <- data[complete.cases(data$Ethnicity_Census),]
data <- data[complete.cases(data$Age_Census),]
data <- data[complete.cases(data$Education_Census),]
data <- data[complete.cases(data$Married_Census),]
data <- data[complete.cases(data$Employment_Status_Census),]

data_C_VP <- data[data$Condition != "Identity_Prime",]
data_C_IP <- data[data$Condition != "Violence_Prime",]
data_IP_VP <- data[data$Condition != "Control",]

data_C_VP$Treatment <- ifelse(data_C_VP$Condition == "Violence_Prime", 1, 0)
data_C_IP$Treatment <- ifelse(data_C_IP$Condition == "Identity_Prime", 1, 0)
data_IP_VP$Treatment <- ifelse(data_IP_VP$Condition == "Violence_Prime", 1, 0)

# CREATE SURVEY WEIGHTS FOR data_C_VP (CONTROL AND VIOLENCE PRIME)
data_C_VP.svy.unweighted <- svydesign(ids=~1, data=data_C_VP)
# Warning message indicates that the initial survey data are not weighted. We create the survey weights manually using the code below.
gender.dist <- data.frame(Male_Census = c(1, 0), Freq = nrow(data_C_VP) * c(0.490566978, 0.509433022))
ethnicity.dist <- data.frame(Ethnicity_Census = c("Bosniak", "Croat", "Serb"), Freq = nrow(data_C_VP) * c(0.520299138, 0.160177354, 0.319523508))
age.dist <- data.frame(Age_Census = c(1, 2, 3, 4), Freq = nrow(data_C_VP) * c(0.202269683, 0.353700847, 0.26717429, 0.17685518))
education.dist <- data.frame(Education_Census = c(1, 2, 3, 4, 5), Freq = nrow(data_C_VP) * c(0.048897384, 0.091729374, 0.214473261, 0.510524395, 0.134375586))
married.dist <- data.frame(Married_Census = c(1, 0), Freq = nrow(data_C_VP) * c(0.588843291, 0.411156709))
employment.dist <- data.frame(Employment_Status_Census = c("Employed", "Unemployed", "Not in labor force"), Freq = nrow(data_C_VP) * c(0.346076909, 0.110004552, 0.543918539))
data_C_VP.svy.rake <- rake(design = data_C_VP.svy.unweighted, 
                        sample.margins = list(~Male_Census, ~Ethnicity_Census, ~Age_Census, ~Education_Census, ~Married_Census, ~Employment_Status_Census), 
                        population.margins = list(gender.dist, ethnicity.dist, age.dist, education.dist, married.dist, employment.dist))
data_C_VP.svy.rake.trim <- trimWeights(data_C_VP.svy.rake, lower=0.3, upper=3, strict=TRUE)
data_C_VP$Weights <- weights(data_C_VP.svy.rake.trim)

# CREATE SURVEY WEIGHTS FOR data_C_IP (CONTROL AND IDENTITY PRIME)
data_C_IP.svy.unweighted <- svydesign(ids=~1, data=data_C_IP)
# Warning message indicates that the initial survey data are not weighted. We create the survey weights manually using the code below.
gender.dist <- data.frame(Male_Census = c(1, 0), Freq = nrow(data_C_IP) * c(0.490566978, 0.509433022))
ethnicity.dist <- data.frame(Ethnicity_Census = c("Bosniak", "Croat", "Serb"), Freq = nrow(data_C_IP) * c(0.520299138, 0.160177354, 0.319523508))
age.dist <- data.frame(Age_Census = c(1, 2, 3, 4), Freq = nrow(data_C_IP) * c(0.202269683, 0.353700847, 0.26717429, 0.17685518))
education.dist <- data.frame(Education_Census = c(1, 2, 3, 4, 5), Freq = nrow(data_C_IP) * c(0.048897384, 0.091729374, 0.214473261, 0.510524395, 0.134375586))
married.dist <- data.frame(Married_Census = c(1, 0), Freq = nrow(data_C_IP) * c(0.588843291, 0.411156709))
employment.dist <- data.frame(Employment_Status_Census = c("Employed", "Unemployed", "Not in labor force"), Freq = nrow(data_C_IP) * c(0.346076909, 0.110004552, 0.543918539))
data_C_IP.svy.rake <- rake(design = data_C_IP.svy.unweighted, 
                        sample.margins = list(~Male_Census, ~Ethnicity_Census, ~Age_Census, ~Education_Census, ~Married_Census, ~Employment_Status_Census), 
                        population.margins = list(gender.dist, ethnicity.dist, age.dist, education.dist, married.dist, employment.dist))
data_C_IP.svy.rake.trim <- trimWeights(data_C_IP.svy.rake, lower=0.3, upper=3, strict=TRUE)
data_C_IP$Weights <- weights(data_C_IP.svy.rake.trim)

# CREATE SURVEY WEIGHTS FOR data_IP_VP (IDENTITY PRIME AND VIOLENCE PRIME)
data_IP_VP.svy.unweighted <- svydesign(ids=~1, data=data_IP_VP)
# Warning message indicates that the initial survey data are not weighted. We create the survey weights manually using the code below.
gender.dist <- data.frame(Male_Census = c(1, 0), Freq = nrow(data_IP_VP) * c(0.490566978, 0.509433022))
ethnicity.dist <- data.frame(Ethnicity_Census = c("Bosniak", "Croat", "Serb"), Freq = nrow(data_IP_VP) * c(0.520299138, 0.160177354, 0.319523508))
age.dist <- data.frame(Age_Census = c(1, 2, 3, 4), Freq = nrow(data_IP_VP) * c(0.202269683, 0.353700847, 0.26717429, 0.17685518))
education.dist <- data.frame(Education_Census = c(1, 2, 3, 4, 5), Freq = nrow(data_IP_VP) * c(0.048897384, 0.091729374, 0.214473261, 0.510524395, 0.134375586))
married.dist <- data.frame(Married_Census = c(1, 0), Freq = nrow(data_IP_VP) * c(0.588843291, 0.411156709))
employment.dist <- data.frame(Employment_Status_Census = c("Employed", "Unemployed", "Not in labor force"), Freq = nrow(data_IP_VP) * c(0.346076909, 0.110004552, 0.543918539))
data_IP_VP.svy.rake <- rake(design = data_IP_VP.svy.unweighted, 
                        sample.margins = list(~Male_Census, ~Ethnicity_Census, ~Age_Census, ~Education_Census, ~Married_Census, ~Employment_Status_Census), 
                        population.margins = list(gender.dist, ethnicity.dist, age.dist, education.dist, married.dist, employment.dist))
data_IP_VP.svy.rake.trim <- trimWeights(data_IP_VP.svy.rake, lower=0.3, upper=3, strict=TRUE)
data_IP_VP$Weights <- weights(data_IP_VP.svy.rake.trim)

# CREATE FINAL SURVEY DESIGN OBJECT
design_C_VP <- svydesign(ids=~1, weights=~Weights, data=data_C_VP)
design_C_IP <- svydesign(ids=~1, weights=~Weights, data=data_C_IP)
design_IP_VP <- svydesign(ids=~1, weights=~Weights, data=data_IP_VP)

# Table SA9.7: Safety Threat and Policy Preferences, OLS Models, Using Survey Weights
# Model 1
summary(SA9.7.1 <- svyglm(Policing_OLS ~ Treatment, design_C_VP))
# extract the number of observations
nobs(SA9.7.1)   
# extract the AIC
SA9.7.1$aic     

# Model 2
summary(SA9.7.2 <- svyglm(Policing_OLS ~ Treatment, design_C_IP))
# extract the number of observations
nobs(SA9.7.2)   
# extract the AIC
SA9.7.2$aic     

# Model 3
summary(SA9.7.3 <- svyglm(Policing_OLS ~ Treatment, design_IP_VP))
# extract the number of observations
nobs(SA9.7.3)   
# extract the AIC
SA9.7.3$aic     

# Table SA9.8: Safety Threat and Policy Preferences, Ordered Probit Models, Using Survey Weights
# Model 1
summary(SA9.8.1 <- svyolr(Policing_OP ~ Treatment, design_C_VP, method = "probit"))
# extract the p-values
coeftest(SA9.8.1)
# extract the number of observations
nobs(SA9.8.1)   
# extract the AIC
summary(SA9.8.1);SA9.8.1$deviance+2*length(SA9.8.1$coefficients) 

# Model 2
summary(SA9.8.2 <- svyolr(Policing_OP ~ Treatment, design_C_IP, method = "probit"))
# extract the p-values
coeftest(SA9.8.2)
# extract the number of observations
nobs(SA9.8.2)   
# extract the AIC
summary(SA9.8.2);SA9.8.2$deviance+2*length(SA9.8.2$coefficients) 

# Model 3
summary(SA9.8.3 <- svyolr(Policing_OP ~ Treatment, design_IP_VP, method = "probit"))
# extract the p-values
coeftest(SA9.8.3)
# extract the number of observations
nobs(SA9.8.3) 
# extract the AIC
summary(SA9.8.3);SA9.8.3$deviance+2*length(SA9.8.3$coefficients) 

# Table SA9.9: Safety Threat and Policy Preferences, Probit Models, Using Survey Weights
# Model 1
summary(SA9.9.1 <- svyglm(Policing_Probit ~ Treatment, design_C_VP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.9.1) 
# extract the AIC
summary(SA9.9.1);SA9.9.1$deviance+2*length(SA9.9.1$coefficients) 

# Model 2
summary(SA9.9.2 <- svyglm(Policing_Probit ~ Treatment, design_C_IP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.9.2) 
# extract the AIC
summary(SA9.9.2);SA9.9.2$deviance+2*length(SA9.9.2$coefficients) 

# Model 3
summary(SA9.9.3 <- svyglm(Policing_Probit ~ Treatment, design_IP_VP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.9.3)  
# extract the AIC
summary(SA9.9.3);SA9.9.3$deviance+2*length(SA9.9.3$coefficients) 

# Table SA9.10: Safety Threat and Policy Preferences, OLS Models, Controlling for Respondent Characteristics, Using Survey Weights
# Model 1
summary(SA9.10.1 <- svyglm(Policing_OLS ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_C_VP))
# extract the number of observations
nobs(SA9.10.1)   
# extract the AIC
SA9.10.1$aic     

# Model 2
summary(SA9.10.2 <- svyglm(Policing_OLS ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_C_VP))
# extract the number of observations
nobs(SA9.10.2)  
# extract the AIC
SA9.10.2$aic     

# Model 3
summary(SA9.10.3 <- svyglm(Policing_OLS ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_C_IP))
# extract the number of observations
nobs(SA9.10.3)   
# extract the AIC
SA9.10.3$aic    

# Model 4
summary(SA9.10.4 <- svyglm(Policing_OLS ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_C_IP))
# extract the number of observations
nobs(SA9.10.4)   
# extract the AIC
SA9.10.4$aic     

# Model 5
summary(SA9.10.5 <- svyglm(Policing_OLS ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_IP_VP))
# extract the number of observations
nobs(SA9.10.5)  
# extract the AIC
SA9.10.5$aic     

# Model 6
summary(SA9.10.6 <- svyglm(Policing_OLS ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_IP_VP))
# extract the number of observations
nobs(SA9.10.6)  
# extract the AIC
SA9.10.6$aic    

# Table SA9.11: Safety Threat and Policy Preferences, Ordered Probit Models, Controlling for Respondent Characteristics, Using Survey Weights
# Model 1
summary(SA9.11.1 <- svyolr(Policing_OP ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_C_VP, method = "probit"))
# extract the p-values
coeftest(SA9.11.1)
# extract the number of observations
nobs(SA9.11.1)   
# extract the AIC
summary(SA9.11.1);SA9.11.1$deviance+2*length(SA9.11.1$coefficients) 

# Model 2
summary(SA9.11.2 <- svyolr(Policing_OP ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_C_VP, method = "probit"))
# extract the p-values
coeftest(SA9.11.2)
# extract the number of observations
nobs(SA9.11.2)  
# extract the AIC
summary(SA9.11.2);SA9.11.2$deviance+2*length(SA9.11.2$coefficients) 

# Model 3
summary(SA9.11.3 <- svyolr(Policing_OP ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_C_IP, method = "probit"))
# extract the p-values
coeftest(SA9.11.3)
# extract the number of observations
nobs(SA9.11.3)  
# extract the AIC
summary(SA9.11.3);SA9.11.3$deviance+2*length(SA9.11.3$coefficients) 

# Model 4
summary(SA9.11.4 <- svyolr(Policing_OP ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_C_IP, method = "probit"))
# extract the p-values
coeftest(SA9.11.4)
# extract the number of observations
nobs(SA9.11.4)   
# extract the AIC
summary(SA9.11.4);SA9.11.4$deviance+2*length(SA9.11.4$coefficients) 

# Model 5
summary(SA9.11.5 <- svyolr(Policing_OP ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_IP_VP, method = "probit"))
# extract the p-values
coeftest(SA9.11.5)
# extract the number of observations
nobs(SA9.11.5)   
# extract the AIC
summary(SA9.11.5);SA9.11.5$deviance+2*length(SA9.11.5$coefficients) 

# Model 6
summary(SA9.11.6 <- svyolr(Policing_OP ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_IP_VP, method = "probit"))
# extract the p-values
coeftest(SA9.11.6)
# extract the number of observations
nobs(SA9.11.6)   
# extract the AIC
summary(SA9.11.6);SA9.11.6$deviance+2*length(SA9.11.6$coefficients) 

# Table SA9.12: Safety Threat and Policy Preferences, Probit Models, Controlling for Respondent Characteristics, Using Survey Weights
# Model 1
summary(SA9.12.1 <- svyglm(Policing_Probit ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_C_VP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.12.1)  
# extract the AIC
summary(SA9.12.1);SA9.12.1$deviance+2*length(SA9.12.1$coefficients) 

# Model 2
summary(SA9.12.2 <- svyglm(Policing_Probit ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_C_VP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.12.2)  
# extract the AIC
summary(SA9.12.2);SA9.12.2$deviance+2*length(SA9.12.2$coefficients)

# Model 3
summary(SA9.12.3 <- svyglm(Policing_Probit ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_C_IP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.12.3)
# extract the AIC
summary(SA9.12.3);SA9.12.3$deviance+2*length(SA9.12.3$coefficients) 

# Model 4
summary(SA9.12.4 <- svyglm(Policing_Probit ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_C_IP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.12.4)  
# extract the AIC
summary(SA9.12.4);SA9.12.4$deviance+2*length(SA9.12.4$coefficients) 

# Model 5
summary(SA9.12.5 <- svyglm(Policing_Probit ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census, design_IP_VP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.12.5)
# extract the AIC
summary(SA9.12.5);SA9.12.5$deviance+2*length(SA9.12.5$coefficients) 

# Model 6
summary(SA9.12.6 <- svyglm(Policing_Probit ~ Treatment + Ethnicity_Census + Male_Census + Age_Census + Education_Census + Married_Census + Employment_Status_Census, design_IP_VP, family=quasibinomial(link="probit")))
# extract the number of observations
nobs(SA9.12.6)  
# extract the AIC
summary(SA9.12.6);SA9.12.6$deviance+2*length(SA9.12.6$coefficients) 


